library(dplyr)
library(raster)
library(rhdf5)
library(rgdal)
library(neonAOP)
library(reshape2)
library(ggplot2)
library(MASS)
library(stargazer)

options(stringsAsFactors = FALSE)

Import D17-SOAP Field Data

## read in file
insitu.SOAP<-read.csv("/Users/R-Files/Documents/data/NEONDI-2016/NEONdata/D17-California/SOAP/2013/insitu/veg-structure/D17_2013_SOAP_vegStr.csv")

## visualize structure
str(insitu.SOAP)
## 'data.frame':    1231 obs. of  30 variables:
##  $ siteid               : chr  "SOAP" "SOAP" "SOAP" "SOAP" ...
##  $ sitename             : chr  "Soaproot Saddle" "Soaproot Saddle" "Soaproot Saddle" "Soaproot Saddle" ...
##  $ plotid               : chr  "SOAP43" "SOAP43" "SOAP43" "SOAP43" ...
##  $ easting              : num  297065 297070 297056 297067 297057 ...
##  $ northing             : num  4100726 4100722 4100722 4100719 4100729 ...
##  $ taxonid              : chr  "ABCO" "ABCO" "ABCO" "ABCO" ...
##  $ scientificname       : chr  "Abies concolor" "Abies concolor" "Abies concolor" "Abies concolor" ...
##  $ indvidual_id         : int  622 632 567 606 639 566 527 638 516 555 ...
##  $ pointid              : chr  "center" "NE" "center" "center" ...
##  $ individualdistance   : num  5.4 5.4 4.1 8.5 5.7 4.1 8.3 6.9 6.4 6 ...
##  $ individualazimuth    : num  67.8 188.7 250.8 122.4 332.6 ...
##  $ dbh                  : num  67.6 16.1 5.9 0.9 4 5.1 1.1 3.1 1 8.9 ...
##  $ dbhheight            : num  130 130 130 10 130 130 130 130 10 130 ...
##  $ basalcanopydiam      : num  8.2 0 0 0 0 0 0 0 0 0 ...
##  $ basalcanopydiam_90deg: num  7.9 0 0 0 0 0 0 0 0 0 ...
##  $ maxcanopydiam        : num  0 3.3 2.9 0.4 2 2.4 1.3 2 0.5 3 ...
##  $ canopydiam_90deg     : num  0 3.2 2.8 0.3 1.5 1.9 1.4 1.9 0.5 2.9 ...
##  $ stemheight           : num  34.5 13 3.8 0.6 3.4 3 1.4 3.1 0.5 5 ...
##  $ stemremarks          : chr  "" "" "" "" ...
##  $ stemstatus           : chr  "" "" "" "" ...
##  $ canopyform           : chr  "" "" "" "" ...
##  $ livingcanopy         : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ inplotcanopy         : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ materialsampleid     : chr  "" "" "" "" ...
##  $ dbhqf                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ stemmapqf            : int  0 1 0 0 0 0 0 0 0 0 ...
##  $ plant_group          : logi  NA NA NA NA NA NA ...
##  $ common_name          : logi  NA NA NA NA NA NA ...
##  $ aop_plot             : logi  NA NA NA NA NA NA ...
##  $ unique_id            : logi  NA NA NA NA NA NA ...
#summarize
summary(insitu.SOAP)
##     siteid            sitename            plotid             easting      
##  Length:1231        Length:1231        Length:1231        Min.   :297034  
##  Class :character   Class :character   Class :character   1st Qu.:297064  
##  Mode  :character   Mode  :character   Mode  :character   Median :298705  
##                                                           Mean   :298489  
##                                                           3rd Qu.:299809  
##                                                           Max.   :300153  
##                                                                           
##     northing         taxonid          scientificname      indvidual_id   
##  Min.   :4100083   Length:1231        Length:1231        Min.   :   1.0  
##  1st Qu.:4100556   Class :character   Class :character   1st Qu.: 304.5  
##  Median :4100844   Mode  :character   Mode  :character   Median : 612.0  
##  Mean   :4100817                                         Mean   : 632.2  
##  3rd Qu.:4101027                                         3rd Qu.: 920.5  
##  Max.   :4101489                                         Max.   :1341.0  
##                                                                          
##    pointid          individualdistance individualazimuth      dbh        
##  Length:1231        Min.   : 0.000     Min.   :  0.00    Min.   :  0.00  
##  Class :character   1st Qu.: 4.700     1st Qu.: 66.95    1st Qu.:  1.20  
##  Mode  :character   Median : 6.900     Median :160.80    Median :  3.65  
##                     Mean   : 6.903     Mean   :173.57    Mean   : 10.10  
##                     3rd Qu.: 9.100     3rd Qu.:286.30    3rd Qu.: 13.68  
##                     Max.   :64.000     Max.   :359.80    Max.   :159.80  
##                                                          NA's   :249     
##    dbhheight     basalcanopydiam   basalcanopydiam_90deg maxcanopydiam   
##  Min.   :  0.0   Min.   :  0.000   Min.   : 0.000        Min.   : 0.000  
##  1st Qu.: 10.0   1st Qu.:  0.000   1st Qu.: 0.000        1st Qu.: 0.700  
##  Median :130.0   Median :  0.000   Median : 0.000        Median : 1.600  
##  Mean   : 86.3   Mean   :  7.556   Mean   : 5.322        Mean   : 2.355  
##  3rd Qu.:130.0   3rd Qu.:  4.000   3rd Qu.: 2.000        3rd Qu.: 3.100  
##  Max.   :140.0   Max.   :220.000   Max.   :93.000        Max.   :80.000  
##                  NA's   :1                               NA's   :1       
##  canopydiam_90deg   stemheight      stemremarks         stemstatus       
##  Min.   : 0.000   Min.   : -1.600   Length:1231        Length:1231       
##  1st Qu.: 0.500   1st Qu.:  1.000   Class :character   Class :character  
##  Median : 1.300   Median :  1.800   Mode  :character   Mode  :character  
##  Mean   : 1.863   Mean   :  4.786                                        
##  3rd Qu.: 2.400   3rd Qu.:  4.450                                        
##  Max.   :70.000   Max.   :134.000                                        
##  NA's   :1                                                               
##   canopyform         livingcanopy     inplotcanopy    materialsampleid  
##  Length:1231        Min.   :  0.00   Min.   : 40.00   Length:1231       
##  Class :character   1st Qu.: 60.00   1st Qu.:100.00   Class :character  
##  Mode  :character   Median :100.00   Median :100.00   Mode  :character  
##                     Mean   : 74.25   Mean   : 99.78                     
##                     3rd Qu.:100.00   3rd Qu.:100.00                     
##                     Max.   :100.00   Max.   :100.00                     
##                     NA's   :85       NA's   :84                         
##      dbhqf     stemmapqf       plant_group    common_name   
##  Min.   :0   Min.   :0.00000   Mode:logical   Mode:logical  
##  1st Qu.:0   1st Qu.:0.00000   NA's:1231      NA's:1231     
##  Median :0   Median :0.00000                                
##  Mean   :0   Mean   :0.09748                                
##  3rd Qu.:0   3rd Qu.:0.00000                                
##  Max.   :0   Max.   :1.00000                                
##  NA's   :1                                                  
##  aop_plot       unique_id     
##  Mode:logical   Mode:logical  
##  NA's:1231      NA's:1231     
##                               
##                               
##                               
##                               
## 

Exploratory Analysis: Subset variables of interest and identify data type (Cat & Cont)

#  selected potential predictors of interest at SOAP 2013

# plotid, taxonid, scientificname, dbh, dbhheight, basalCanopyDiam, Maximum canopy diameter, stemheight
# canopyform, livingcanopy,  inplotcanopy

# subset dataframe-- make it not dirty
insitu.SOAP<-subset(insitu.SOAP, select=c(plotid, easting, northing, taxonid, scientificname, dbh, dbhheight, basalcanopydiam, maxcanopydiam, stemheight, canopyform, livingcanopy,  inplotcanopy))

## subset categorical variables
SOAP.cat<-subset(insitu.SOAP, select=c(plotid, taxonid, scientificname, canopyform))

## subset continuous variables
SOAP.cont<- subset(insitu.SOAP, select=c(dbh, dbhheight, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy))

Exploratory Analysis for SOAP- Frequency Categorical Statistics

## frequency summaries by:

# plotid
SOAP.cat %>%
  group_by(plotid) %>%
  summarise (n = n()) %>%
  mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))
## Source: local data frame [18 x 3]
## 
##      plotid     n rel.freq
##       (chr) (int)    (chr)
## 1  SOAP1343    24     1.9%
## 2   SOAP139   122     9.9%
## 3   SOAP143   114     9.3%
## 4  SOAP1515    26     2.1%
## 5  SOAP1563    14     1.1%
## 6  SOAP1611    10     0.8%
## 7  SOAP1695   106     8.6%
## 8   SOAP187    95     7.7%
## 9   SOAP223    43     3.5%
## 10  SOAP255   104     8.4%
## 11  SOAP283    14     1.1%
## 12  SOAP299    26     2.1%
## 13  SOAP331   132    10.7%
## 14   SOAP43   100     8.1%
## 15  SOAP555     6     0.5%
## 16   SOAP63   122     9.9%
## 17   SOAP95   106     8.6%
## 18  SOAP991    67     5.4%
# dominant(?) plant type

SOAP.cat %>%
  group_by(scientificname) %>%
  summarise (n = n()) %>%
  mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))
## Source: local data frame [16 x 3]
## 
##                      scientificname     n rel.freq
##                               (chr) (int)    (chr)
## 1                    Abies concolor    11     0.9%
## 2             Amelanchier utahensis     2     0.2%
## 3            Arctostaphylos viscida   179    14.5%
## 4              Calocedrus decurrens   435    35.3%
## 5                Ceanothus cuneatus    49       4%
## 6            Ceanothus integerrimus    42     3.4%
## 7  Cercocarpus montanus var. glaber    40     3.2%
## 8                   Corylus cornuta     3     0.2%
## 9                 Pinus lambertiana    16     1.3%
## 10                  Pinus ponderosa   213    17.3%
## 11                Prunus emarginata     1     0.1%
## 12              Quercus chrysolepis   141    11.5%
## 13                 Quercus garryana     3     0.2%
## 14                Quercus kelloggii    65     5.3%
## 15               Rhamnus ilicifolia    12       1%
## 16                    Ribes roezlii    19     1.5%
# canopyform

SOAP.cat %>%
  group_by(canopyform) %>%
  summarise (n = n()) %>%
  mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))
## Source: local data frame [6 x 3]
## 
##   canopyform     n rel.freq
##        (chr) (int)    (chr)
## 1              844    68.6%
## 2       Cone   155    12.6%
## 3   Cylinder    14     1.1%
## 4 Hemisphere    70     5.7%
## 5     Sphere    64     5.2%
## 6         NA    84     6.8%

Exploratory Analysis for SOAP- Summary Continuous Statistics

## Continuous variable summaries- Min/Max, Quants
summary(SOAP.cont)
##       dbh           dbhheight     basalcanopydiam   maxcanopydiam   
##  Min.   :  0.00   Min.   :  0.0   Min.   :  0.000   Min.   : 0.000  
##  1st Qu.:  1.20   1st Qu.: 10.0   1st Qu.:  0.000   1st Qu.: 0.700  
##  Median :  3.65   Median :130.0   Median :  0.000   Median : 1.600  
##  Mean   : 10.10   Mean   : 86.3   Mean   :  7.556   Mean   : 2.355  
##  3rd Qu.: 13.68   3rd Qu.:130.0   3rd Qu.:  4.000   3rd Qu.: 3.100  
##  Max.   :159.80   Max.   :140.0   Max.   :220.000   Max.   :80.000  
##  NA's   :249                      NA's   :1         NA's   :1       
##    stemheight       livingcanopy     inplotcanopy   
##  Min.   : -1.600   Min.   :  0.00   Min.   : 40.00  
##  1st Qu.:  1.000   1st Qu.: 60.00   1st Qu.:100.00  
##  Median :  1.800   Median :100.00   Median :100.00  
##  Mean   :  4.786   Mean   : 74.25   Mean   : 99.78  
##  3rd Qu.:  4.450   3rd Qu.:100.00   3rd Qu.:100.00  
##  Max.   :134.000   Max.   :100.00   Max.   :100.00  
##                    NA's   :85       NA's   :84
## Continuous variable summaries- Std Devs
sapply(SOAP.cont, sd, na.rm=TRUE)
##             dbh       dbhheight basalcanopydiam   maxcanopydiam 
##       14.927547       58.329815       19.304295        3.803459 
##      stemheight    livingcanopy    inplotcanopy 
##        8.526687       41.145058        2.863307
# matrix scatterplot
plot(SOAP.cont)

#side by side boxplots
par(mar=c(9,5,1,1))
boxplot(SOAP.cont, las=2)

## boxplot by station ID
dat<-subset(x = insitu.SOAP, select= c(plotid, dbh, basalcanopydiam, maxcanopydiam, stemheight), na.rm=TRUE)

par(mfrow=c(2,2))
boxplot(dbh~plotid, data=dat, las=2, main="dbh (cm)")
boxplot(basalcanopydiam~plotid, data=dat, las=2, main="basalcanopydiam (cm)")
boxplot(maxcanopydiam~plotid, data=dat, las=2, main="maxcanopydiam (m)")
boxplot(stemheight~plotid, data=dat, las=2, main="stemheight (m)")

Exploratory Analysis for SOAP by Sampling Plots

## by species types by PlotID

xtabs(~plotid+scientificname, data=SOAP.cat)
##           scientificname
## plotid     Abies concolor Amelanchier utahensis Arctostaphylos viscida
##   SOAP1343              0                     0                     12
##   SOAP139               0                     0                     14
##   SOAP143               0                     0                     24
##   SOAP1515              0                     0                     23
##   SOAP1563              0                     0                     14
##   SOAP1611              0                     0                     10
##   SOAP1695              0                     0                     36
##   SOAP187               0                     0                      0
##   SOAP223               0                     0                      8
##   SOAP255               0                     0                     15
##   SOAP283               0                     0                      7
##   SOAP299               0                     0                      3
##   SOAP331               4                     2                      0
##   SOAP43                7                     0                      0
##   SOAP555               0                     0                      6
##   SOAP63                0                     0                      6
##   SOAP95                0                     0                      0
##   SOAP991               0                     0                      1
##           scientificname
## plotid     Calocedrus decurrens Ceanothus cuneatus Ceanothus integerrimus
##   SOAP1343                    1                  5                      0
##   SOAP139                    72                  0                      1
##   SOAP143                    13                  0                     34
##   SOAP1515                    0                  0                      1
##   SOAP1563                    0                  0                      0
##   SOAP1611                    0                  0                      0
##   SOAP1695                    0                  1                      0
##   SOAP187                    81                  0                      0
##   SOAP223                     0                  2                      0
##   SOAP255                    14                  0                      0
##   SOAP283                     3                  0                      1
##   SOAP299                     3                  0                      0
##   SOAP331                    85                  0                      0
##   SOAP43                     63                  0                      0
##   SOAP555                     0                  0                      0
##   SOAP63                     38                  0                      0
##   SOAP95                     61                  0                      5
##   SOAP991                     1                 41                      0
##           scientificname
## plotid     Cercocarpus montanus var. glaber Corylus cornuta
##   SOAP1343                                2               0
##   SOAP139                                 0               0
##   SOAP143                                 0               0
##   SOAP1515                                0               0
##   SOAP1563                                0               0
##   SOAP1611                                0               0
##   SOAP1695                                0               0
##   SOAP187                                 1               0
##   SOAP223                                16               0
##   SOAP255                                 0               0
##   SOAP283                                 0               0
##   SOAP299                                 0               0
##   SOAP331                                 0               0
##   SOAP43                                  0               3
##   SOAP555                                 0               0
##   SOAP63                                  0               0
##   SOAP95                                  0               0
##   SOAP991                                21               0
##           scientificname
## plotid     Pinus lambertiana Pinus ponderosa Prunus emarginata
##   SOAP1343                 0               1                 0
##   SOAP139                  3               6                 0
##   SOAP143                  0              38                 0
##   SOAP1515                 0               1                 0
##   SOAP1563                 0               0                 0
##   SOAP1611                 0               0                 0
##   SOAP1695                 0               4                 0
##   SOAP187                  0               8                 0
##   SOAP223                  0               0                 0
##   SOAP255                  0              51                 0
##   SOAP283                  0               0                 0
##   SOAP299                  0              16                 0
##   SOAP331                  1              11                 0
##   SOAP43                   9               6                 0
##   SOAP555                  0               0                 0
##   SOAP63                   0              53                 0
##   SOAP95                   3              18                 0
##   SOAP991                  0               0                 1
##           scientificname
## plotid     Quercus chrysolepis Quercus garryana Quercus kelloggii
##   SOAP1343                   0                1                 1
##   SOAP139                   19                0                 2
##   SOAP143                    0                0                 1
##   SOAP1515                   1                0                 0
##   SOAP1563                   0                0                 0
##   SOAP1611                   0                0                 0
##   SOAP1695                  50                0                 0
##   SOAP187                    2                0                 2
##   SOAP223                   14                1                 1
##   SOAP255                   15                0                 9
##   SOAP283                    3                0                 0
##   SOAP299                    4                0                 0
##   SOAP331                    6                0                23
##   SOAP43                     0                0                 9
##   SOAP555                    0                0                 0
##   SOAP63                    10                0                15
##   SOAP95                    17                0                 2
##   SOAP991                    0                1                 0
##           scientificname
## plotid     Rhamnus ilicifolia Ribes roezlii
##   SOAP1343                  0             1
##   SOAP139                   0             5
##   SOAP143                   0             4
##   SOAP1515                  0             0
##   SOAP1563                  0             0
##   SOAP1611                  0             0
##   SOAP1695                 12             3
##   SOAP187                   0             1
##   SOAP223                   0             1
##   SOAP255                   0             0
##   SOAP283                   0             0
##   SOAP299                   0             0
##   SOAP331                   0             0
##   SOAP43                    0             3
##   SOAP555                   0             0
##   SOAP63                    0             0
##   SOAP95                    0             0
##   SOAP991                   0             1
## aggregate summary stats 
#  Re-instate PlotID

SOAP.cont$plotid<-insitu.SOAP$plotid

# group by
melted <- melt(SOAP.cont, id.vars=c("plotid"))

SOAP.aggregate<-summarise(group_by(melted, plotid, variable),
          mean=mean(value, na.rm=TRUE), sd=sd(value, na.rm=TRUE), 
          min=min(value, na.rm=TRUE), max=max(value, na.rm=TRUE))
 
SOAP.aggregate
## Source: local data frame [126 x 6]
## Groups: plotid [?]
## 
##      plotid        variable       mean        sd   min   max
##       (chr)          (fctr)      (dbl)     (dbl) (dbl) (dbl)
## 1  SOAP1343             dbh   8.400000  6.437391   1.0  15.4
## 2  SOAP1343       dbhheight 120.000000 33.879582  10.0 130.0
## 3  SOAP1343 basalcanopydiam  29.287500 30.840844   0.0 140.0
## 4  SOAP1343   maxcanopydiam   3.920833  5.736684   0.6  30.0
## 5  SOAP1343      stemheight   2.112500  1.212368   0.7   6.5
## 6  SOAP1343    livingcanopy  71.625000 33.217089   0.0 100.0
## 7  SOAP1343    inplotcanopy 100.000000  0.000000 100.0 100.0
## 8   SOAP139             dbh   7.274510 19.207526   0.4 159.8
## 9   SOAP139       dbhheight  80.983607 59.089648  10.0 130.0
## 10  SOAP139 basalcanopydiam   5.311475 13.882674   0.0  60.0
## ..      ...             ...        ...       ...   ...   ...
## WHOA-- That's a big spicy meatball!
## Break it down by plot again-- this time stratify by parameter

#overside columnn width restriction
options(dplyr.width = Inf)

# dbh
 insitu_dbh <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(dbh.min = min(dbh, na.rm=TRUE), dbh.max = max(dbh, na.rm=TRUE), dbh.mean=mean(dbh, na.rm=TRUE), dbh.median=median(dbh, na.rm=TRUE),dbh.sd=sd(dbh, na.rm=TRUE))

 insitu_dbh             
## Source: local data frame [18 x 6]
## 
##      plotid dbh.min dbh.max  dbh.mean dbh.median    dbh.sd
##       (chr)   (dbl)   (dbl)     (dbl)      (dbl)     (dbl)
## 1  SOAP1343     1.0    15.4  8.400000       8.60  6.437391
## 2   SOAP139     0.4   159.8  7.274510       1.70 19.207526
## 3   SOAP143     0.4    44.0  6.018889       1.40 11.503276
## 4  SOAP1515     1.0    24.0 11.646154      10.50  6.357753
## 5  SOAP1563     0.0     0.0  0.000000       0.00  0.000000
## 6  SOAP1611     0.0     0.0  0.000000       0.00  0.000000
## 7  SOAP1695     0.0    60.7  7.654667       1.90 12.417656
## 8   SOAP187     1.2   102.9  8.220652       2.75 16.507677
## 9   SOAP223     0.6   107.9 16.823529      10.20 26.028387
## 10  SOAP255     0.4    78.3 10.751724       7.30 13.063271
## 11  SOAP283     0.0    33.5  6.492857       0.00 12.273578
## 12  SOAP299     1.1    64.1 32.685000      33.90 21.016767
## 13  SOAP331     0.4    72.1 10.714729       5.90 12.610971
## 14   SOAP43     0.7    84.0 11.066316       5.80 16.032312
## 15  SOAP555     0.0     0.0  0.000000       0.00  0.000000
## 16   SOAP63     0.2    50.8 14.746364      13.25 12.035816
## 17   SOAP95     0.3    61.3 10.762626       5.50 12.969474
## 18  SOAP991     1.4    10.5  6.420000       8.20  4.490212
## functionalization might be best here...but for now 
 ## JV scripting

 
  # dbhheight
 insitu_dbhheight <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(dbhheight.min = min(dbhheight, na.rm=TRUE), dbhheight.max = max(dbhheight, na.rm=TRUE), dbhheight.mean=mean(dbhheight, na.rm=TRUE), dbhheight.median=median(dbhheight, na.rm=TRUE),dbhheight.sd=sd(dbhheight, na.rm=TRUE))
 
 insitu_dbhheight
## Source: local data frame [18 x 6]
## 
##      plotid dbhheight.min dbhheight.max dbhheight.mean dbhheight.median
##       (chr)         (dbl)         (dbl)          (dbl)            (dbl)
## 1  SOAP1343          10.0           130      120.00000            130.0
## 2   SOAP139          10.0           130       80.98361            130.0
## 3   SOAP143          10.0           130       59.47368             10.0
## 4  SOAP1515           1.0           130       70.88462            130.0
## 5  SOAP1563           0.0             0        0.00000              0.0
## 6  SOAP1611           0.0             0        0.00000              0.0
## 7  SOAP1695           0.0           130       66.36604             38.5
## 8   SOAP187           2.0           130       69.07368             10.0
## 9   SOAP223          10.0           130       94.18605            130.0
## 10  SOAP255          10.0           140      109.32692            130.0
## 11  SOAP283           0.0           130       28.18571              0.0
## 12  SOAP299          10.0           130      125.38462            130.0
## 13  SOAP331          10.0           130       83.65909            130.0
## 14   SOAP43          10.0           130      103.60000            130.0
## 15  SOAP555           0.0             0        0.00000              0.0
## 16   SOAP63           0.8           130      111.23607            130.0
## 17   SOAP95          10.0           130       85.84906            130.0
## 18  SOAP991          10.0           130      121.04478            130.0
##    dbhheight.sd
##           (dbl)
## 1      33.87958
## 2      59.08965
## 3      59.33022
## 4      65.13944
## 5       0.00000
## 6       0.00000
## 7      61.59260
## 8      60.22498
## 9      55.12896
## 10     45.67221
## 11     55.19133
## 12     23.53394
## 13     58.62428
## 14     49.95998
## 15      0.00000
## 16     43.87607
## 17     58.14379
## 18     31.77260
 # basalcanopydiam
 insitu_basalcanopydiam <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(basalcanopydiam.min = min(basalcanopydiam, na.rm=TRUE), basalcanopydiam.max = max(basalcanopydiam, na.rm=TRUE), basalcanopydiam.mean=mean(basalcanopydiam, na.rm=TRUE), basalcanopydiam.median=median(basalcanopydiam, na.rm=TRUE),basalcanopydiam.sd=sd(basalcanopydiam, na.rm=TRUE))
 
 insitu_basalcanopydiam
## Source: local data frame [18 x 6]
## 
##      plotid basalcanopydiam.min basalcanopydiam.max basalcanopydiam.mean
##       (chr)               (dbl)               (dbl)                (dbl)
## 1  SOAP1343                 0.0                 140           29.2875000
## 2   SOAP139                 0.0                  60            5.3114754
## 3   SOAP143                 0.0                  40            2.3771930
## 4  SOAP1515                 0.0                 117           23.8461538
## 5  SOAP1563                 2.0                  90           27.2142857
## 6  SOAP1611                 1.1                 220           42.6100000
## 7  SOAP1695                 0.0                  94           10.3000000
## 8   SOAP187                 0.0                  48            1.0315789
## 9   SOAP223                 0.0                  58           11.3255814
## 10  SOAP255                 0.0                  80            5.5000000
## 11  SOAP283                 0.0                  31           14.0000000
## 12  SOAP299                 0.0                 125           15.9230769
## 13  SOAP331                 0.0                  60            0.5174242
## 14   SOAP43                 0.0                  38            0.8320000
## 15  SOAP555                 7.0                 130           90.8333333
## 16   SOAP63                 0.0                  61            2.3360656
## 17   SOAP95                 0.0                  30            0.7809524
## 18  SOAP991                 0.0                 163           34.6567164
##    basalcanopydiam.median basalcanopydiam.sd
##                     (dbl)              (dbl)
## 1                    22.5          30.840844
## 2                     0.0          13.882674
## 3                     0.0           6.012362
## 4                    20.0          21.502916
## 5                    24.0          22.874839
## 6                    18.0          64.730286
## 7                     0.0          17.087255
## 8                     0.0           6.104502
## 9                     6.0          15.085030
## 10                    0.0          15.270379
## 11                   13.0          11.948608
## 12                    0.0          35.305436
## 13                    0.0           5.253389
## 14                    0.0           4.755664
## 15                  110.0          48.803347
## 16                    0.0           8.759987
## 17                    0.0           3.855518
## 18                   24.0          31.471754
 # maxcanopydiam
 insitu_maxcanopydiam <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(maxcanopydiam.min = min(maxcanopydiam, na.rm=TRUE), maxcanopydiam.max = max(maxcanopydiam, na.rm=TRUE), maxcanopydiam.mean=mean(maxcanopydiam, na.rm=TRUE), maxcanopydiam.median=median(maxcanopydiam, na.rm=TRUE),maxcanopydiam.sd=sd(maxcanopydiam, na.rm=TRUE))
 
 insitu_maxcanopydiam
## Source: local data frame [18 x 6]
## 
##      plotid maxcanopydiam.min maxcanopydiam.max maxcanopydiam.mean
##       (chr)             (dbl)             (dbl)              (dbl)
## 1  SOAP1343               0.6              30.0           3.920833
## 2   SOAP139               0.2              30.0           1.752459
## 3   SOAP143               0.2              11.5           1.240351
## 4  SOAP1515               0.7               5.4           2.334615
## 5  SOAP1563               0.5               3.9           1.885714
## 6  SOAP1611               1.0               5.1           2.240000
## 7  SOAP1695               0.1              60.0           3.513208
## 8   SOAP187               0.0              80.0           2.309474
## 9   SOAP223               0.3              10.1           2.576744
## 10  SOAP255               0.0              12.6           3.031731
## 11  SOAP283               0.0               3.8           1.807143
## 12  SOAP299               0.0              10.1           4.330769
## 13  SOAP331               0.0               8.5           1.934351
## 14   SOAP43               0.0              18.1           2.440000
## 15  SOAP555               1.3               4.0           2.200000
## 16   SOAP63               0.0               6.8           2.604918
## 17   SOAP95               0.0              12.4           2.096226
## 18  SOAP991               0.4               4.7           1.962687
##    maxcanopydiam.median maxcanopydiam.sd
##                   (dbl)            (dbl)
## 1                  2.85        5.7366843
## 2                  0.90        3.2085684
## 3                  0.60        1.8008709
## 4                  2.10        1.0403624
## 5                  1.80        1.0021954
## 6                  1.80        1.3150412
## 7                  2.20        7.1522634
## 8                  0.80        8.3154186
## 9                  1.90        2.1445209
## 10                 2.40        2.3197449
## 11                 1.70        1.2256642
## 12                 4.25        2.9714332
## 13                 1.40        1.6596928
## 14                 1.95        2.4516743
## 15                 1.95        0.9423375
## 16                 2.45        1.6878688
## 17                 0.90        2.2582304
## 18                 2.00        0.9766727
  # stemheight
 insitu_stemheight <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(stemheight.min = min(stemheight, na.rm=TRUE), stemheight.max = max(stemheight, na.rm=TRUE), stemheight.mean=mean(stemheight, na.rm=TRUE), stemheight.median=median(stemheight, na.rm=TRUE),stemheight.sd=sd(stemheight, na.rm=TRUE))
 
 insitu_stemheight
## Source: local data frame [18 x 6]
## 
##      plotid stemheight.min stemheight.max stemheight.mean
##       (chr)          (dbl)          (dbl)           (dbl)
## 1  SOAP1343            0.7            6.5       2.1125000
## 2   SOAP139            0.5          120.0       4.7459016
## 3   SOAP143            0.5           19.7       2.5307018
## 4  SOAP1515            0.6            3.4       1.7769231
## 5  SOAP1563            0.0            1.9       1.0714286
## 6  SOAP1611            0.6            4.2       1.7200000
## 7  SOAP1695            0.6           17.0       3.0566038
## 8   SOAP187            0.5          134.0       4.8789474
## 9   SOAP223            0.5            8.2       2.4441860
## 10  SOAP255            0.0           44.7       5.1019231
## 11  SOAP283            0.0           13.8       3.3285714
## 12  SOAP299            1.0           33.1      13.8961538
## 13  SOAP331            0.5           39.2       5.0098485
## 14   SOAP43            0.5           51.1       5.0470000
## 15  SOAP555            0.8            1.4       0.9666667
## 16   SOAP63            0.6           33.0       9.2327869
## 17   SOAP95            0.5           28.1       6.0509434
## 18  SOAP991           -1.6            3.8       1.8537313
##    stemheight.median stemheight.sd
##                (dbl)         (dbl)
## 1               2.05     1.2123683
## 2               1.50    13.0803565
## 3               1.00     4.1863057
## 4               1.70     0.5867249
## 5               1.20     0.5426836
## 6               1.30     1.2318008
## 7               2.20     2.8116770
## 8               1.40    14.9537862
## 9               1.60     1.9816182
## 10              2.35     6.7463328
## 11              1.80     4.1795157
## 12             10.30    12.8098550
## 13              2.40     6.4031226
## 14              2.95     7.8090010
## 15              0.90     0.2250926
## 16              5.00     9.0170403
## 17              1.70     8.2020583
## 18              1.90     0.9241422
  # livingcanopy
 insitu_livingcanopy <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(livingcanopy.min = min(livingcanopy, na.rm=TRUE), livingcanopy.max = max(livingcanopy, na.rm=TRUE), livingcanopy.mean=mean(livingcanopy, na.rm=TRUE), livingcanopy.median=median(livingcanopy, na.rm=TRUE),livingcanopy.sd=sd(livingcanopy, na.rm=TRUE))
  
 insitu_livingcanopy
## Source: local data frame [18 x 6]
## 
##      plotid livingcanopy.min livingcanopy.max livingcanopy.mean
##       (chr)            (int)            (int)             (dbl)
## 1  SOAP1343                0              100          71.62500
## 2   SOAP139                0              100          79.71311
## 3   SOAP143                0              100          95.92105
## 4  SOAP1515               20              100          82.42308
## 5  SOAP1563               NA               NA               NaN
## 6  SOAP1611               NA               NA               NaN
## 7  SOAP1695                0              100          88.57576
## 8   SOAP187                0              100          23.15789
## 9   SOAP223                0              100          79.48837
## 10  SOAP255                0              100          66.89320
## 11  SOAP283               NA               NA               NaN
## 12  SOAP299                0              100          84.61538
## 13  SOAP331                0              100          78.78788
## 14   SOAP43                0              100          81.90000
## 15  SOAP555               NA               NA               NaN
## 16   SOAP63                0              100          62.58197
## 17   SOAP95                0              100          84.90566
## 18  SOAP991                0              100          71.41791
##    livingcanopy.median livingcanopy.sd
##                  (dbl)           (dbl)
## 1                   80        33.21709
## 2                  100        39.42405
## 3                  100        18.47703
## 4                   90        20.29418
## 5                   NA             NaN
## 6                   NA             NaN
## 7                  100        19.37172
## 8                    0        42.40793
## 9                   90        24.98225
## 10                 100        46.67953
## 11                  NA             NaN
## 12                 100        36.79465
## 13                 100        41.03676
## 14                 100        38.57814
## 15                  NA             NaN
## 16                 100        48.22568
## 17                 100        35.96944
## 18                  80        26.73957
   # inplotcanopy
 insitu_inplotcanopy <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(inplotcanopy.min = min(inplotcanopy, na.rm=TRUE), inplotcanopy.max = max(inplotcanopy, na.rm=TRUE), inplotcanopy.mean=mean(inplotcanopy, na.rm=TRUE), inplotcanopy.median=median(inplotcanopy, na.rm=TRUE),inplotcanopy.sd=sd(inplotcanopy, na.rm=TRUE))
 
 insitu_inplotcanopy 
## Source: local data frame [18 x 6]
## 
##      plotid inplotcanopy.min inplotcanopy.max inplotcanopy.mean
##       (chr)            (int)            (int)             (dbl)
## 1  SOAP1343              100              100         100.00000
## 2   SOAP139               90              100          99.91803
## 3   SOAP143              100              100         100.00000
## 4  SOAP1515              100              100         100.00000
## 5  SOAP1563               NA               NA               NaN
## 6  SOAP1611               NA               NA               NaN
## 7  SOAP1695               65              100          98.03030
## 8   SOAP187              100              100         100.00000
## 9   SOAP223              100              100         100.00000
## 10  SOAP255              100              100         100.00000
## 11  SOAP283               NA               NA               NaN
## 12  SOAP299              100              100         100.00000
## 13  SOAP331              100              100         100.00000
## 14   SOAP43               40              100          98.90000
## 15  SOAP555               NA               NA               NaN
## 16   SOAP63              100              100         100.00000
## 17   SOAP95              100              100         100.00000
## 18  SOAP991              100              100         100.00000
##    inplotcanopy.median inplotcanopy.sd
##                  (dbl)           (dbl)
## 1                  100       0.0000000
## 2                  100       0.9053575
## 3                  100       0.0000000
## 4                  100       0.0000000
## 5                   NA             NaN
## 6                   NA             NaN
## 7                  100       6.7867965
## 8                  100       0.0000000
## 9                  100       0.0000000
## 10                 100       0.0000000
## 11                  NA             NaN
## 12                 100       0.0000000
## 13                 100       0.0000000
## 14                 100       7.7713538
## 15                  NA             NaN
## 16                 100       0.0000000
## 17                 100       0.0000000
## 18                 100       0.0000000

Let’s Extract Some Vegetation Indices for SOAP

## import NDVI
SOAP.ndvi <- raster("../NEONdata/D17-California/SOAP/2013/spectrometer/veg_index/SOAP_NDVI.tif")
SOAP.ndvi
## class       : RasterLayer 
## dimensions  : 1516, 3292, 4990672  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 296906, 300198, 4100038, 4101554  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/R-Files/Documents/data/NEONDI-2016/NEONdata/D17-California/SOAP/2013/spectrometer/veg_index/SOAP_NDVI.tif 
## names       : SOAP_NDVI 
## values      : -0.5165877, 1  (min, max)
##visualize distribution

hist(SOAP.ndvi, main="NDVI for SOAP Site")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 2% of the raster cells were used. 100000 values used.

# visualize raster
plot(SOAP.ndvi, main="NDVI for SOAP Site")

## read in plot centroids
SOAP_plots <- readOGR("../NEONdata/D17-California/SOAP/vector_data/",
                      "SOAP_centroids")
## OGR data source with driver: ESRI Shapefile 
## Source: "../NEONdata/D17-California/SOAP/vector_data/", layer: "SOAP_centroids"
## with 50 features
## It has 10 fields
# Overlay the centroid points and the stem locations on the NDVI plot
plot(SOAP.ndvi,
     main="SOAP Plot Locations \n(N=18)",
     col=gray.colors(100, start=.3, end=.9))

# pch 0 = square
plot(SOAP_plots,
     pch = 0,
     cex = 2,
     col = 4,
     add=TRUE)

# Insitu sampling took place within 40m x 40m square plots, so we use a 20m radius.
# Note that below will return a dataframe containing the max height
# calculated from all pixels in the buffer for each plot (NEON Website)

## we will use these values to populate in-situ data set

NDVI.max <- extract(SOAP.ndvi,
                    SOAP_plots,
                    buffer = 20,
                    fun=(max),
                    sp=TRUE,
                    stringsAsFactors=FALSE)


# NDVI distribution visualization for each plot
# created histograms for all 14 plots of data

cent_ovrList <- extract(SOAP.ndvi,SOAP_plots,buffer = 20)

for (i in 1:14) {
  hist(cent_ovrList[[i]], main=(paste("plot",i)))
  }

Extract descriptive stats from Insitu Data

# import the centroid data and the vegetation structure data

## Call out insitu.SOAP dataframe
## refamilize with PlotIDs within SOAP
unique(insitu.SOAP$plotid)
##  [1] "SOAP43"   "SOAP331"  "SOAP139"  "SOAP1343" "SOAP143"  "SOAP63"  
##  [7] "SOAP1563" "SOAP1695" "SOAP255"  "SOAP1611" "SOAP283"  "SOAP1515"
## [13] "SOAP223"  "SOAP555"  "SOAP299"  "SOAP991"  "SOAP95"   "SOAP187"
# find the max parameter values for each in situ plot
insitu.max <- insitu.SOAP %>%
  group_by(plotid) %>%
  summarise(max.dbh = max(dbh, na.rm=TRUE), max.dbhheight=max(dbhheight, na.rm=TRUE),
            max.basalcanopydiam= max(basalcanopydiam, na.rm=TRUE), 
            max.maxcanopydiam=max(maxcanopydiam), 
            max.stemheight=max(stemheight, na.rm=TRUE),
            max.livingcanopy=max(livingcanopy, na.rm=TRUE),
            max.inplotcanopy= max(inplotcanopy, na.rm=TRUE))
## merge data into data frame, can use approach to create spatial object
## But purposes of capstone will conduct "boring regression analysis"
## merge grouped data (insitu.max) and extracted NDVI (NDVI.max) 
## using plotid to merge

## first must append  NDVI.max file with SOAP

#create data frame for NDVI.max

NDVI.max.df<-NDVI.max@data

# add SOAP prefix to IDs

NDVI.max.df$plot.id<- "SOAP"
NDVI.max.df$plot.id<-(sapply(NDVI.max.df$plot.id, paste0, NDVI.max.df$ID))[,1]

# rename column

colnames(NDVI.max.df)[colnames(NDVI.max.df) == 'SOAP'] <- 'Plot.ID'

SOAP.grand<-merge(NDVI.max.df, insitu.max, 
      by.x = 'plot.id',
      by.y = 'plotid')

# F* Yes
## Reponse variables: dbh, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy

## Predictor variable: SOAP_NDVI

## Simple Linear Regression

## Max NDVI vs Max DBH
mod.ndvi<-lm(max.dbh ~ SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.dbh, SOAP.grand$SOAP_NDVI, xlab="Max NDVI", ylab=" DBH (cm)")
abline(mod.ndvi)

## Max NDVI vs Max Base Canopy Diam
mod2.ndvi<-lm(max.basalcanopydiam~ SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.basalcanopydiam, SOAP.grand$SOAP_NDVI, 
     xlab="Max NDVI", ylab=" Base Canopy Diameter (cm)")
abline(mod2.ndvi)

## Max NDVI vs Max of the Max Canopy Diam
mod3.ndvi<-lm(max.maxcanopydiam~SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.maxcanopydiam, SOAP.grand$SOAP_NDVI, 
     xlab="Max NDVI", ylab=" Maximum canopy diameter(m)")
abline(mod3.ndvi)

## Max NDVI vs Max Stem Height
mod4.ndvi<-lm(max.stemheight~SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.stemheight, SOAP.grand$SOAP_NDVI, 
     xlab="Max NDVI", ylab=" Stem Height (m)")
abline(mod4.ndvi)

create a qualitative landscape assessment index (completely experimental- BE NICE!)

#subset SOAP.grand data set

#s<-subset(SOAP.grand, select = c(plot.id, easting, northing, max.dbh, max.dbhheight, max.basalcanopydiam, max.maxcanopydiam, max.stemheight))

## categorize by quartiles, use quantile values from summary and assign score

# dbh
summary(SOAP.grand$max.dbh)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   41.38   62.70   65.91   88.72  159.80
dbh.cat <- cut(SOAP.grand$max.dbh, breaks = c(-1, 41.38, 62.70, 88.72, 159.8))
dbh.cat<-unclass(dbh.cat)


#basalcanopydiam
summary(SOAP.grand$max.basalcanopydiam)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   30.00   39.50   59.00   77.83   77.00  220.00
max.basalcanopydiam.cat<- cut(SOAP.grand$max.basalcanopydiam, breaks= c(0, 39.50, 59.0, 77.00, 222.00))
max.basalcanopydiam.cat<- unclass(max.basalcanopydiam.cat)

# maxcanopydiam
summary(SOAP.grand$max.maxcanopydiam)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    3.80    5.95   10.10   17.51   15.25   80.00       1
## impute missing value at Plot SOAP331
SOAP.grand$max.maxcanopydiam[8] <- mean(SOAP.grand$max.maxcanopydiam, na.rm=TRUE)

# determine class quantiles scores
summary(SOAP.grand$max.maxcanopydiam)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.800   6.375  10.800  17.510  17.660  80.000
max.maxcanopydiam.cat<-cut(SOAP.grand$max.maxcanopydiam, breaks= c(0, 6.375, 10.80, 17.660, 80.00))
max.maxcanopydiam.cat<-unclass(max.maxcanopydiam.cat)

#stem height
summary(SOAP.grand$max.stemheight)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3.80   12.40   30.55   40.68   42.18  134.00
max.stemheight.cat<-cut(SOAP.grand$max.stemheight, breaks= c(0,12.40,30.55, 42.18, 134))
max.stemheight.cat<-unclass(max.stemheight.cat)

s<-data.frame(dbh.cat, max.basalcanopydiam.cat,max.maxcanopydiam.cat,max.stemheight.cat) 

colnames(s)<-c("dbh", "basalcanopydial", "maxcanopydiam", "stemheight")

tree.score<-rowSums(s)

s$tree.score<-tree.score
s$tree.score.perc<- (tree.score/16)

n<-data.frame(c(SOAP.grand$plot.id),s)

colnames(n)= c("plotid","dbh", "basalcanopydial", "maxcanopydiam", "stemheight", "Score", "% Score")

## polished score table
stargazer(n, type="html", out="score.table.htm")
## 
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">dbh</td><td>12</td><td>2.500</td><td>1.168</td><td>1</td><td>4</td></tr>
## <tr><td style="text-align:left">basalcanopydial</td><td>12</td><td>2.500</td><td>1.168</td><td>1</td><td>4</td></tr>
## <tr><td style="text-align:left">maxcanopydiam</td><td>12</td><td>2.500</td><td>1.168</td><td>1</td><td>4</td></tr>
## <tr><td style="text-align:left">stemheight</td><td>12</td><td>2.500</td><td>1.168</td><td>1</td><td>4</td></tr>
## <tr><td style="text-align:left">Score</td><td>12</td><td>10.000</td><td>3.045</td><td>5</td><td>15</td></tr>
## <tr><td style="text-align:left">% Score</td><td>12</td><td>0.625</td><td>0.190</td><td>0.312</td><td>0.938</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>
stargazer(n, type="text", out="score.table.txt")
## 
## ==============================================
## Statistic       N   Mean  St. Dev.  Min   Max 
## ----------------------------------------------
## dbh             12 2.500   1.168     1     4  
## basalcanopydial 12 2.500   1.168     1     4  
## maxcanopydiam   12 2.500   1.168     1     4  
## stemheight      12 2.500   1.168     1     4  
## Score           12 10.000  3.045     5    15  
## % Score         12 0.625   0.190   0.312 0.938
## ----------------------------------------------
## add on to SOAP.grand dataset

SOAP.grand<-c(SOAP.grand,s)

boxplot(s$tree.score.perc)

## simple linear regression using tree score ()

mod5.ndvi<-lm(tree.score.perc~SOAP_NDVI, data=SOAP.grand)

## also- makesure you append CHM values from other data stack to model
## import chm
SOAP.chm <- raster("../NEONdata/D17-California/SOAP/2013/lidar/SOAP_lidarCHM.tif")
SOAP.chm[SOAP.chm==0]<- NA
plot(SOAP.chm, main="Canopy Height Model for SOAP Site")

# Insitu sampling took place within 40m x 40m square plots, so we use a 20m radius.
# Note that below will return a dataframe containing the max height
# calculated from all pixels in the buffer for each plot (NEON Website)

## we will use these values to populate in-situ data set

chm.max <- extract(SOAP.chm,
                    SOAP_plots,
                    buffer = 20,
                    fun=(max),
                    sp=TRUE,
                    stringsAsFactors=FALSE)


# CHM distribution visualization for each plot
# created histograms for all 14 plots of data

cent_ovrList <- extract(SOAP.chm,SOAP_plots,buffer = 20)

for (i in 1:14) {
  hist(cent_ovrList[[i]], main=(paste("plot",i)))
  }

## merge data into data frame, can use approach to create spatial object
## But purposes of capstone will conduct "boring regression analysis"
## merge grouped data (insitu.max) and extracted NDVI (NDVI.max) 
## using plotid to merge

## first must append  NDVI.max file with SOAP

#create data frame for NDVI.max

chm.max.df<-chm.max@data

# add SOAP prefix to IDs

chm.max.df$plot.id<- "SOAP"
chm.max.df$plot.id<-(sapply(chm.max.df$plot.id, paste0, chm.max.df$ID))[,1]

# rename column

colnames(chm.max.df)[colnames(chm.max.df) == 'SOAP'] <- 'Plot.ID'

SOAP.grand<-merge(chm.max.df, SOAP.grand, 
      by.x = 'plot.id',
      by.y = 'plot.id')

# F* Yes
## Response variables: dbh, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy

## Predictor variable: SOAP_lidarCHM

## Simple Linear Fit

mod.chm<-lm(max.dbh~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.dbh, SOAP.grand$SOAP_lidarCHM, 
     xlab="Max CHM", ylab=" DBH (cm)")
abline(mod.chm)

summary(mod.chm)
## 
## Call:
## lm(formula = max.dbh ~ SOAP_lidarCHM, data = SOAP.grand)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.517 -16.469  -8.646   0.459  78.467 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     4.9555    30.4248   0.163    0.874  
## SOAP_lidarCHM   1.8172     0.8449   2.151    0.057 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.35 on 10 degrees of freedom
## Multiple R-squared:  0.3163, Adjusted R-squared:  0.2479 
## F-statistic: 4.626 on 1 and 10 DF,  p-value: 0.05698
mod2.chm<-lm(max.basalcanopydiam~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.basalcanopydiam, SOAP.grand$SOAP_lidarCHM, 
     xlab="Max CHM", ylab=" Base Canopy Diameter (cm)")
abline(mod2.chm)

summary(mod2.chm)
## 
## Call:
## lm(formula = max.basalcanopydiam ~ SOAP_lidarCHM, data = SOAP.grand)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.543 -47.516  -5.578  19.236 110.410 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    133.947     46.005   2.912   0.0155 *
## SOAP_lidarCHM   -1.673      1.278  -1.309   0.2197  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 57.99 on 10 degrees of freedom
## Multiple R-squared:  0.1464, Adjusted R-squared:  0.06102 
## F-statistic: 1.715 on 1 and 10 DF,  p-value: 0.2197
mod3.chm<-lm(max.maxcanopydiam~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.maxcanopydiam, SOAP.grand$SOAP_lidarCHM, 
     xlab="Max CHM", ylab=" Maximum canopy diameter(m)")
abline(mod3.chm)

summary(mod3.chm)
## 
## Call:
## lm(formula = max.maxcanopydiam ~ SOAP_lidarCHM, data = SOAP.grand)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.240  -9.587  -3.554   0.542  49.986 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -11.5842    14.4213  -0.803   0.4405  
## SOAP_lidarCHM   0.8674     0.4005   2.166   0.0556 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.18 on 10 degrees of freedom
## Multiple R-squared:  0.3193, Adjusted R-squared:  0.2512 
## F-statistic: 4.691 on 1 and 10 DF,  p-value: 0.05556
mod4.chm<-lm(max.stemheight~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.stemheight, SOAP.grand$SOAP_lidarCHM, 
     xlab="Max CHM", ylab=" Stem Height (m)")
abline(mod4.chm)

summary(mod4.chm)
## 
## Call:
## lm(formula = max.stemheight ~ SOAP_lidarCHM, data = SOAP.grand)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.832 -17.010  -5.808   9.709  59.479 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -38.0415    23.7283  -1.603  0.13997   
## SOAP_lidarCHM   2.3470     0.6589   3.562  0.00516 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 29.91 on 10 degrees of freedom
## Multiple R-squared:  0.5592, Adjusted R-squared:  0.5151 
## F-statistic: 12.69 on 1 and 10 DF,  p-value: 0.005165
## add on to SOAP.grand dataset

SOAP.grand<-c(SOAP.grand,s)
boxplot(s$tree.score.perc)

## simple linear regression using tree score 

mod5.chm<-lm(tree.score.perc~SOAP_lidarCHM, data=SOAP.grand)
mod5.chm
## 
## Call:
## lm(formula = tree.score.perc ~ SOAP_lidarCHM, data = SOAP.grand)
## 
## Coefficients:
##   (Intercept)  SOAP_lidarCHM  
##       0.24781        0.01125
## Full model in respect to dbh
mod.full<- lm(max.dbh~SOAP_lidarCHM + SOAP_NDVI, data=SOAP.grand)
mod.full
## 
## Call:
## lm(formula = max.dbh ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
## 
## Coefficients:
##   (Intercept)  SOAP_lidarCHM      SOAP_NDVI  
##      -137.057          1.505        172.771
## Full model in respect to base canopy diameter
mod2.full<-lm(max.basalcanopydiam~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod2.full
## 
## Call:
## lm(formula = max.basalcanopydiam ~ SOAP_lidarCHM + SOAP_NDVI, 
##     data = SOAP.grand)
## 
## Coefficients:
##   (Intercept)  SOAP_lidarCHM      SOAP_NDVI  
##     856.50684       -0.08263     -879.05838
## Full model in respect to max canopy
mod3.full<-lm(max.maxcanopydiam~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod3.full
## 
## Call:
## lm(formula = max.maxcanopydiam ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
## 
## Coefficients:
##   (Intercept)  SOAP_lidarCHM      SOAP_NDVI  
##      -15.1373         0.8595         4.3227
## Full model in respect to stem height
mod4.full<-lm(max.stemheight~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod4.full
## 
## Call:
## lm(formula = max.stemheight ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
## 
## Coefficients:
##   (Intercept)  SOAP_lidarCHM      SOAP_NDVI  
##        -9.646          2.410        -34.546
## Full model in respect to % score
mod5.full<- lm(tree.score.perc~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod5.full
## 
## Call:
## lm(formula = tree.score.perc ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
## 
## Coefficients:
##   (Intercept)  SOAP_lidarCHM      SOAP_NDVI  
##       0.41638        0.01162       -0.20508
#descriptive statistics for in-situ field data products
stargazer(SOAP.cont, type = "text", title="SOAP-wide Descriptive statistics", digits=1, out="SOAPwideStats.txt")
## 
## SOAP-wide Descriptive statistics
## ==============================================
## Statistic         N   Mean St. Dev. Min   Max 
## ----------------------------------------------
## dbh              982  10.1   14.9   0.0  159.8
## dbhheight       1,231 86.3   58.3   0.0  140.0
## basalcanopydiam 1,230 7.6    19.3   0.0  220.0
## maxcanopydiam   1,230 2.4    3.8    0.0  80.0 
## stemheight      1,231 4.8    8.5    -1.6 134.0
## livingcanopy    1,146 74.2   41.1    0    100 
## inplotcanopy    1,147 99.8   2.9     40   100 
## ----------------------------------------------
stargazer(data.frame(insitu.max), type = "text", title="Maximum Value Plot-wide Descriptive statistics", digits=1, out="PlotwideStats-freq.txt")
## 
## Maximum Value Plot-wide Descriptive statistics
## ===============================================
## Statistic           N  Mean  St. Dev. Min  Max 
## -----------------------------------------------
## max.dbh             18 53.9    43.4   0.0 159.8
## max.dbhheight       18 108.9   50.2    0   140 
## max.basalcanopydiam 18 88.1    52.0   30   220 
## max.maxcanopydiam   17 18.1    21.4   3.8 80.0 
## max.stemheight      18 31.3    38.3   1.4 134.0
## max.livingcanopy    14 100.0   0.0    100  100 
## max.inplotcanopy    14 100.0   0.0    100  100 
## -----------------------------------------------
## regression products

## NDVI-only

stargazer(mod.ndvi, mod2.ndvi, mod3.ndvi, mod4.ndvi, mod5.ndvi, type="html", dep.var.labels=c("dbh,cm", "basalCanopyDiam,cm", "maxCanopyDiam,m", "stemHeight,m", "Combinatory Tree Index(-)"),covariate.labels=("NDVI"), ci=TRUE, note= "NDVI as the only predictor. Maximum NDVI values determined from 20 ft buffers at each respective SOAP station ", out="NDVImodels2.htm")
## 
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="5"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="5" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>dbh,cm</td><td>basalCanopyDiam,cm</td><td>maxCanopyDiam,m</td><td>stemHeight,m</td><td>Combinatory Tree Index(-)</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">NDVI</td><td>394.748</td><td>-891.248<sup>***</sup></td><td>132.966</td><td>320.930</td><td>1.509</td></tr>
## <tr><td style="text-align:left"></td><td>(-122.167, 911.663)</td><td>(-1,433.297, -349.200)</td><td>(-142.008, 407.940)</td><td>(-197.733, 839.592)</td><td>(-0.768, 3.785)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-282.516</td><td>864.495<sup>***</sup></td><td>-100.070</td><td>-242.586</td><td>-0.707</td></tr>
## <tr><td style="text-align:left"></td><td>(-739.388, 174.355)</td><td>(385.409, 1,343.581)</td><td>(-343.575, 143.435)</td><td>(-701.002, 215.831)</td><td>(-2.719, 1.306)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>12</td><td>12</td><td>11</td><td>12</td><td>12</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.183</td><td>0.509</td><td>0.091</td><td>0.128</td><td>0.144</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.101</td><td>0.460</td><td>-0.010</td><td>0.041</td><td>0.059</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>41.920 (df = 10)</td><td>43.959 (df = 10)</td><td>22.145 (df = 9)</td><td>42.062 (df = 10)</td><td>0.185 (df = 10)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>2.240 (df = 1; 10)</td><td>10.385<sup>***</sup> (df = 1; 10)</td><td>0.898 (df = 1; 9)</td><td>1.471 (df = 1; 10)</td><td>1.687 (df = 1; 10)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="5" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
## 
## <table style="text-align:center"><tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr><tr><td>NDVI as the only predictor. Maximum NDVI values determined from 20 ft buffers at each respective SOAP station</td></tr>
## <tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr></table>
## CHM-only
stargazer(mod.chm, mod2.chm, mod3.chm, mod4.chm, mod5.chm, type="html", dep.var.labels=c("dbh,cm", "basalCanopyDiam,cm", "maxCanopyDiam,m", "stemHeight,m", "Combinatory Tree Index(-)"),covariate.labels=("CHM"), ci=TRUE, note= "CHM as the only predictor. Maximum CHM values determined from 20 ft buffers at each respective SOAP station ", out="CHMmodels2.htm")
## 
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="5"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="5" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>dbh,cm</td><td>basalCanopyDiam,cm</td><td>maxCanopyDiam,m</td><td>stemHeight,m</td><td>Combinatory Tree Index(-)</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">CHM</td><td>1.817<sup>*</sup></td><td>-1.673</td><td>0.867<sup>*</sup></td><td>2.347<sup>***</sup></td><td>0.011<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.161, 3.473)</td><td>(-4.177, 0.831)</td><td>(0.082, 1.652)</td><td>(1.056, 3.638)</td><td>(0.006, 0.016)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>4.955</td><td>133.947<sup>**</sup></td><td>-11.584</td><td>-38.042</td><td>0.248<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(-54.676, 64.587)</td><td>(43.779, 224.116)</td><td>(-39.849, 16.681)</td><td>(-84.548, 8.465)</td><td>(0.065, 0.430)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>12</td><td>12</td><td>12</td><td>12</td><td>12</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.316</td><td>0.146</td><td>0.319</td><td>0.559</td><td>0.654</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.248</td><td>0.061</td><td>0.251</td><td>0.515</td><td>0.619</td></tr>
## <tr><td style="text-align:left">Residual Std. Error (df = 10)</td><td>38.349</td><td>57.988</td><td>18.177</td><td>29.909</td><td>0.117</td></tr>
## <tr><td style="text-align:left">F Statistic (df = 1; 10)</td><td>4.626<sup>*</sup></td><td>1.715</td><td>4.691<sup>*</sup></td><td>12.687<sup>***</sup></td><td>18.891<sup>***</sup></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="5" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
## 
## <table style="text-align:center"><tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr><tr><td>CHM as the only predictor. Maximum CHM values determined from 20 ft buffers at each respective SOAP station</td></tr>
## <tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr></table>
## Full (NDVI + CHM)

stargazer(mod.full, mod2.full, mod3.full, mod4.full, mod5.full, type="html", dep.var.labels=c("dbh,cm", "basalCanopyDiam,cm", "maxCanopyDiam,m", "stemHeight,m", "Combinatory Tree Index(-)"),covariate.labels=("NDVI+ CHM"), ci=TRUE, note= "NDVI and CHM as predictors. Maximum NDVI and CHM values determined from 20 ft buffers at each respective SOAP station ", out="Fullmodels.htm2")
## 
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="5"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="5" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>dbh,cm</td><td>basalCanopyDiam,cm</td><td>maxCanopyDiam,m</td><td>stemHeight,m</td><td>Combinatory Tree Index(-)</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">NDVI+ CHM</td><td>1.505</td><td>-0.083</td><td>0.860</td><td>2.410<sup>**</sup></td><td>0.012<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(-0.495, 3.505)</td><td>(-2.419, 2.254)</td><td>(-0.107, 1.826)</td><td>(0.822, 3.997)</td><td>(0.005, 0.018)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">SOAP_NDVI</td><td>172.771</td><td>-879.058<sup>**</sup></td><td>4.323</td><td>-34.546</td><td>-0.205</td></tr>
## <tr><td style="text-align:left"></td><td>(-398.347, 743.888)</td><td>(-1,546.200, -211.917)</td><td>(-271.608, 280.253)</td><td>(-488.015, 418.923)</td><td>(-1.983, 1.573)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-137.057</td><td>856.507<sup>**</sup></td><td>-15.137</td><td>-9.646</td><td>0.416</td></tr>
## <tr><td style="text-align:left"></td><td>(-610.531, 336.417)</td><td>(303.426, 1,409.588)</td><td>(-243.892, 213.618)</td><td>(-385.586, 366.294)</td><td>(-1.057, 1.890)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>12</td><td>12</td><td>12</td><td>12</td><td>12</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.342</td><td>0.510</td><td>0.319</td><td>0.560</td><td>0.656</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.196</td><td>0.401</td><td>0.168</td><td>0.463</td><td>0.579</td></tr>
## <tr><td style="text-align:left">Residual Std. Error (df = 9)</td><td>39.657</td><td>46.324</td><td>19.160</td><td>31.487</td><td>0.123</td></tr>
## <tr><td style="text-align:left">F Statistic (df = 2; 9)</td><td>2.339</td><td>4.678<sup>**</sup></td><td>2.112</td><td>5.735<sup>**</sup></td><td>8.575<sup>***</sup></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="5" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
## 
## <table style="text-align:center"><tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr><tr><td>NDVI and CHM as predictors. Maximum NDVI and CHM values determined from 20 ft buffers at each respective SOAP station</td></tr>
## <tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr></table>